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To estimate cosmological parameters from a given dataset, we need to construct a likelihood func- 
tion, which sometimes has a complicated functional form. We introduce the copula, a mathematical 
tool to construct an arbitrary multivariate distribution function from one-dimensional marginal dis- 
tribution functions with any given dependence structure. It is shown that a likelihood function 
constructed by the so-called Gaussian copula can reproduce very well the n-dimensional probability 
distribution of the cosmic shear power spectrum obtained from a large number of ray-tracing sim- 
ulations. This suggests that the Copula likelihood will be a powerful tool for future weak lensing 
analyses, instead of the conventional multivariate Gaussian likelihood. 

PACS numbers: 98.80.Es 

I. INTRODUCTION 

Weak gravitational lensing by intervening large scale cosmic structures provides an excellent tool to probe the 
nature of dark matter and dark energy. The so-called "cosmic shear" signal has been successfully measured in various 
groups since 2000 [e.g. !TI-fTol| ■ If systematic errors are well under control, the weak lensing has the highest potential 
to constrain the physical parameters in the equation of state governing the dark energy among the cosmological 
observations, such as type la supernovae, baryon acoustic oscillations, and number count of galaxy clusters [Ill4l3| . 

A number of wide-field weak lensing surveys have been planned for this purpose, such as the Subaru Hyper Suprime- 
Cam Survey [l4j], the Panoramic Survey Telescope & Rapid Response System (Pan-STARRS 1 ), the Dark Energy Sur- 
vey (DES 2 ), the Large Synoptic Survey Telescope (LSST 3 ), the Joint Dark Energy Mission (JDEM 4 ) and Euclid 
It is expected that such wide-field weak lensing surveys will reduce statistical errors significantly compared to the past 
and ongoing surveys, because the number of observed galaxies increases proportional to the survey area. However, to 
make a maximal use of the full potential of planned weak lensing surveys for estimating cosmological parameters, it is 
of great importance to employ adequate statistical measures and methods for weak lensing. Particularly, one needs to 
take into account properly the correlations of the observables between different angular scales and redshifts, i.e. the 
covariances. Furthermore, we also need to use an appropriate likelihood function with given marginal distributions. 
If we do not use proper statistical measures and methods, or if we adopt an inaccurate covariance and/or a likelihood 
function, obtained results may be systematically biased 0, [TtJ • 

For the cosmological parameter estimation, almost all previous authors used the x 2 method in weak lensing analyses 
[e.g. [oT-fToj. However, it is found that the probability distribution function (PDF) of the weak lensing power spectrum 
is well approximated by the % 2 distribution though it has a larger positive tail than expected from the % 2 distribution 
[l8j . The x 2 distribution deviates from the Gaussian distribution on large scales because the number of modes 
corresponding to the degree of freedom is very small. Meanwhile, the x 2 distribution approaches the Gaussian 
distribution at high I due to the central limit theorem. We have to include this information accurately when we place 
constraints on the cosmological parameters. 

If all of the marginal distributions are Gaussian distributed, it is straightforward to reconstruct the multivariate 
PDF, which is the so-called multivariate Gaussian PDF. However, it is not a trivial task to reconstruct the original 
multivariate PDF from general marginal distributions. There has been an infinite number of the degree of freedom in 
choosing the original PDF unless the dependence structure is specified. The copula provides us with a straightforward 
solution to this problem. The copula has been used in the field of mathematical finance, but not been widely used in 
the field of astronomy and cosmology, except only a few applications [l9l - [23| . Hence, a copula may have a potential 
to open new fields of astronomy and cosmology. 
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In this paper, we construct a more plausible likelihood function using the Gaussian copula (hereafter "Copula 
likelihood") than the multivariate Gaussian likelihood for the cosmic shear power spectrum. We show that the 
Copula likelihood well reproduces the n-dimensional probability distribution of the cosmic she ar p ower spectrum 
estimated from 1000 realizations which is obtained from ray-tracing simulations performed by |18| . Cosmological 
parameters employed for our ray-tracing simulations are consistent with the WMAP 3-year results (24|. The detail 
descriptions of our ray-tracing simulations are summarized in [l8| (see also, [25[ ). In a companion paper [26j], we 
estimate the cosmological parameters using both the copula likelihood constructed by this paper and the Gaussian 
likelihood in order to evaluate how the difference between two likelihoods affects the parameter estimation. 



II. FORMULATION 



The likelihood function is a central tool for any kind of parameter estimation. It is defined as a function of 
parameters in a given statistical model. The likelihood function £ is related to the joint probability density function 
(JPDF) denoted by f(xi,x 2 , . . . ,x n \9) as 

C(9 | xi,x 2 , . . . , x n ) = f(xi,x 2 , ...,£„ | 0), (1) 

where ii (i = 1, 2, . . . , n) are independent and identically-distributed observed variables and 9 denotes model param- 
eters. We will suppress the argument 9 hereafter. Recall that the n-point cumulative distribution function (CDF), 
denoted by F, is defined as 



F(xi,x 2 ,...,x n ) = / / ... / f(xi,x 2 ,...,x n )dx 1 dx 2 ...dx n , (2) 

J— oo J — oo J —oo 

and there is a following relation between JPDF and CDF: 

* - x d n F(x 1 ,x 2 ,-.-,x n ) 

J{Xi,X 2 , ■ ■ ■ ,X n ) = TT1 • (3) 

ax\Ox 2 ■ ■ ■ ox n 

From the Sklar's theorem [27}, we can obtain the following relation: 

Prob(xi < xi,x 2 <x 2 ,-..,x n <x n ) = F(£i, x 2 , . . . , x n ) = C{F\{x\), F 2 (x 2 ), . . . , F n (x n )), (4) 
where C denotes the function called copula and Fi denotes one-point CDF defined by 

Fi(xi)= / fi{x)dx = Ui . (5) 



Therefore, the copula indicates how the one-point CDFs are jointed together to give the n-point CDF. A comprehensive 
proof of Sklar's theorem and rigorous definition of a copula are found in (23l . 28 1. From the above relation, we can 
easily derive that ii — F[~ (ui), and then derive the following relation, 

C(u 1 ,u 2 ,...,u n )=F(F- 1 (u 1 ),F- 1 (u 2 ),... : F- 1 (u n )) , (6) 

from Eq. (j4j . Differentiating Eq. ([6]) with Eq. ([3]) gives the density of copula c as 

d n C(u 1 ,u 2 ,...,u n ) f(xi,x 2 , . . . ,x n ) 
c{u u u 2 ,...,u n ) = — — = , 7) 

n 

where each fi is the marginal density function of the marginal CDF Fi . The JPDF can then be expressed as 

n 

f(xi,x 2 ,...,x n ) =c(ui,u 2 ,...,u n ) JJ/i(xj). (8) 

1=1 

If the variables Xi are independent of each other, c = 1. In general, however, they are often correlated and c shows 
their correlations as a function of one-point CDF of each stochastic variable. 
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III. GAUSSIAN COPULA 

In this section, we derive the Copula likelihood using a Gaussian copula which is more plausible than the multivariate 
Gaussian likelihood for the cosmic shear power spectrum. The multivariate Gaussian copula is the copula of the n- 
dimensional random vector that is normally distributed. This copula is expressed by 

c( ul , U2 ,...,u n ) = $($^ l («i).*r 1 M.-".*r 1 («n)) • (9) 

Here one-point Gaussian CDF $1 is 



—00 



— ex ( (X ~ M ' )2 



$i(sO=/ -5=- exp ( - ^ n ™> ) dx = ut, (10) 



and $ is n-point Gaussian CDF defined by 



Xl pX 2 



$>(xi,X2,...,x n ) = I I ... I 1 ==cxp ( --(x - /x) T Cov r (x- fi)) dx 1 dx 2 ...dx n) (11) 

, v /(27r)"det(Cov) \ ^ 



00 J —oo 



where we consider the n-point Gaussian CDF with mean /x and n x n covariance matrix. Cov^ 1 shows the inverse 
covariance matrix. We define /x = (pi, [i 2 , ■ ■ ■ , ^n) and x = (x\, x 2 , ■ ■ ■ , x n ) and superscript 'T' stands for the transpose 
of vector. 



A. Gaussian PDF 



Let us consider a two- variables case for simplicity. First consider the case that one-point PDF is Gaussian with <n 
standard deviation and fii mean. In this case, the Gaussian copula is 



C(x u x 2 ) = $ $^ 1 ($ 1 (f 2 ))) , 

1 



(12) 



$(fi,4 2 ) = / / 

J —OO J — 



— OO y/ (27r) 2 det(Cov) 



exp ^— — {x — /x) T Cov x (x — /x)^ dxidx 2 = ^ ^ <fi(xi,x 2 )dxidx 2 . 



From Eq. (J7J, we derive the density of copula as 

9$(ii,x 2 ) 9 2 $ dxidx 2 



c(ui,u 2 ) 



By using 



du\du 2 
duj 



1 



dx\dx 2 du\ du 2 

(it - m) 2 



{x\,x 2 ) 



dx\ dx 2 
du\ du 2 



2na, 



exp 



2a 2 



we finally obtain the density of copula as 

1 



c(ui,u 2 ) 



^/(27r) 2 dot(Cov) 



exp (— ~(x — /x) T Cov 1 (x — //)) 



/27TCT 

0-1(72 



exp 



: exp 



-i^-MrCCov- 1 -^ 3 )- 1 )^-/*)] 



x /det(Cov) 

where J stands for the identity matrix. Therefore, we can obtain the JPDF as 



f(xi,x 2 ) = 



0\0 2 



v /det(Cov) 
1 



exp 



-^[(x-^CCov- 1 -^ 2 )- 1 )^-/^)] 



n 



2na., 



exp 



2a, 2 



exp f (x — zx) T Cov X (A — u) 

v /(2^) 2 det(Cov) \ 2 l p ' v ^; 



(13) 



(14) 



(15) 



(16) 



(17) 



from Eq. ©. Therefore, the JPDF using the Gaussian copula with a Gaussian one-point PDF results in the multi- 
variate Gaussian distribution, as expected. 
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B. Beyond Gaussian PDF 

Now, let us consider the case in which a one-point PDF is not a Gaussian distribution but a general probability 
distribution /j. In this case, the Gaussian copula is 

C{x 1} x 2 ) = $ (^7 1 (F 1 (x 1 )),^ 1 (F 2 (x 2 ))) , (18) 

where one-point CDF Fi is 



F, 



(xi) = / fi(x)dx = Ui. (19) 

J — oo 



Defining as qi = Q 1 1 (Fi(xi)) — ^ 1 1 (ui), we obtain the copula density from Eq. (|18|) as 

c f Ul U2 ) = g2$ (gi'g2) = d 2 3> d qi dq 2 
' du\du 2 dqidq 2 du\ du 2 ' 



where 



dui dui \ dq t ) \\/2iia l \ 2a, 

We use a formula of the differential of the inverse function at the second equality of the above equation. By using 
Eq. (HD) and Eq. flST}, we rewrite Eq. flU as 



m T \u t) _ f^MY 1 _ { i .^/te-^,, . (21) 



-'(tb^K-^))" (vS5-(- 1s ^ £ )) 



From Eq. ((8|), we can easily calculate the JPDF as 



1 ( 1/ ..,Tp„.-l,„ .A/ 1 / (9i -Mi) 



f(xi,X2)= 7== / ==exp --(g-/x) Cov (g - m) 7= exp 



v /(2 7 r) 2 det(Cov) V 2 V V^i V 2cr 

■ — exp ~ tT 



-1 



cxp ( - ^f_™ ) ) h(xi)f 2 (x 2 ). (23) 



We extend the above equation to the n-dimensional case and then finally obtain the n-dimensional JPDF as 
• • • , «n) = V(2ff)n ^ (Cov) «P " MrOov-^g - „)) ft (7^ «P (-^1^) ) 

(24) 

If /i is Gaussian, qi becomes ii. In this case, we can see that f(xi,x 2 , . . . ,x n ) (Eq. [24)) reduces to a multivariate 
Gaussian distribution. Therefore, when fi are not Gaussian PDFs, Eq. ([M]) carries the whole information on the 
correction to the Gaussian distribution. 

Both in practice and from theoretical point of view (29j . it is more appropriate to work with the logarithm of 
the likelihood function, ln£, called log-likelihood. For example, the log-likelihood is used for the likelihood ratio 
test. The likelihood ratio test statistic is defined as twice the difference in these log-likelihoods with a minus sign. 
This quantity is also fundamental for the information statistics and information criterion theory (see, e.g. (29| and 
references therein). From Eq. (|24|) . this test statistic is derived as 

- - - (ft - Hi) 2 " 
-21n£ c (£i,i 2 , ...,x n )= /]/Xqi - ^Cov -1 ^ - fij) - ^ — — 2^ 2 ^ln (25) 

i=l j=l i—1 1 i=l 

for a general probability distribution. For a Gaussian case it reduces to the well-known form 

n n 

- 21n£ g (xi,i 2 , ...,£„) = ^^(^i - //i)Cov _1 (xj - (26) 

i=i j=i 
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Here we abbreviate the irrelevant constant term in the above two equations. 

Finally, we give the relation between &i and qi because one has to calculate given Xj. Since Xi is related to it, 
through Eq. (fT9|) . we have only to derive the relation between Ui and qi. From qi = $^~ 1 (tti), we get 

Ui = = /" -^exp (- {X ~ f i i)2 ) dx. (27) 

We change the variable x to y — (x — /ii)/<7j and then obtain 

' 1 exp(-4)dy. (28) 



2^ 7-oo V 2 



At a first glance, we can recognize this equation as m = \&i [ — — — ) , where Vfi is the cumulative standard normal 
distribution. We can obtain the relation between and % as 

^—^ = ^(ui) % = (T^VO + Mi- (29) 



IV. RESULT AND DISCUSSION 



In this section, we investigate whether the Copula likelihood (Eq. [25]) reproduces the n-dimensional JPDF obtained 
from ray-tracing simulations performed by [18j . We then compare the shape of the Copula likelihood with that of 
the Gaussian likelihood and how the Copula likelihood is statistically better than the Gaussian likelihood. In our 
application, the observed variables X4 are the binned convergence power spectra P K (ii), which are estimated from 
each realization. We estimate it for an assumed bin width Ahai. Throughout this paper, we employ the bin width 
A\n£ = 0.3 and assume that the single source redshift distribution, i.e. all lensed galaxies lie at z s — 1.0 and we do 
not consider intrinsic ellipticity dispersion a e . In our ray-tracing simulations, the survey area is set as fl s = 25 deg 2 . 
Therefore the fundamental mode of our ray-tracing simulations is at the multipole if — 72. We take 13 bins for our 
likelihood analysis. Therefore, we cover the multipole range from if = 72 to i max = 2635. We call the convergence 
power spectrum estimated at bin i as "bin i power" . 

The left panel of Fig.UJshows the two-dimensional JPDF between bin 1 power and bin 2 power from 1000 realizations 
of ray-tracing simulations. The one-point PDF is normalized so that the mean power spectrum is equal to unity. The 
red solid and blue solid contours show the la and 2a confidence level (CL) regions, respectively. In the right panel 
of Fig. [TJ the blue and red contours show two-dimensional marginalized la and 2a confidence regions on the bin 1 
power and bin 2 power plane, which are derived from the Gaussian likelihood (Eq. 121)1) and Copula likelihood (Eq. |2"5|) , 
respectively. One can clearly see that the Copula likelihood model reproduces the simulation data (left panel of Fig. [1]) 
much better than the Gaussian likelihood model, as discussed below. 

In our Copula likelihood model, we have chosen a % 2 distribution for a general one-point probability distribution 
fi. This is a nice choice because the one-point PDF of convergence power spectrum is fairly well described by a x 2 
distribution where the mean and variance are P K (ii) = (P K (£i)) and a 2 (ii) = (P K {ii) 2 ) — P K (£i) 2 , respectively [see 
[T^ |. From [3(|, this % 2 distribution is shown as 

MP^,)) = 4^-(T m) j , (30, 

for P K (ii) > and f x 2 — for P K {ii) < 0. Here, r(x) is the gamma function and we define T = P K (ti) 2 ja 2 (ii) which 
corresponds to the number of independent modes. The covariance matrix is also estimated from the simulation. 

In order to obtain these likelihood contours, we employed the Markov Chain Monte Carlo method [3l|. Assuming 
flat priors for bin i powers, we explored bin i power estimations in the multidimensional space (i.e. 13 dimensional 
space in this case). Eight parallel chains are computed and the convergence test is made based on the Gelman and 
Rubin statistics called "R — 1" statistics 32]. Each of our chains typically has 400,000 points and R — 1 < 0.05 in 
both models. 

We can see that the results between our Copula likelihood function (red contours) and the Gaussian likelihood (blue 
contours) are very different. The difference mainly comes from the one-point PDF, which is taken as a % 2 distribution 
or a Gaussian distribution. The % 2 distribution denoted by Eq. (|30[) deviates from the Gaussian distribution on large 
scales, such as a case considered in Fig. [TJ because the number of modes corresponding to the degrees of freedom is 
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FIG. 1: Left panel: Two-dimensional JPDF between convergence power spectrum estimated at bin 1 and bin 2 among 1000 
realizations obtained from ray-tracing simulations. The bin 1 and bin 2 correspond to multipole I — 72 and £ — 97, respectively. 
The one-point PDF is normalized so that the mean convergence power spectrum estimated at each bin gives unity. The red 
solid and blue solid contours show la and 2a CL, respectively. Right panel: Two-dimensional marginalized constraints on 
convergence power spectrum estimated at bin 1 and bin 2. The red and blue contours show the marginalized constraints (la 
and 2a CL) obtained by Eq. (|25[) based on the Gaussian copula model and Gaussian likelihood of Eq. (|26p , respectively. 




very small. We can also see that our likelihood function well reproduces the results from ray-tracing simulations than 
the Gaussian likelihood. In particular, the results of our Copula likelihood function capture the feature that the value 
which takes the maximum probability deviates from the mean value. The values at which the Copula likelihood takes 
its maximum are (0.619,0.241) at bin 1 and bin 2. 

Figure [5] is same as Fig. [TJ but results from bin 6 power and bin 7 power. The bin 6 and bin 7 correspond to 
the multipoles I = 323 and £ = 436, respectively. The results from our likelihood are similar to the results from the 
Gaussian likelihood, because the \ 2 distribution approaches Gaussian distribution at these scales. However, we can 
see the small deviation from each mean value between the results obtained from the two likelihoods. 

Figure [3] is also the same as Fig. [TJ but results from bin 12 power and bin 13 power. The bin 12 and bin 13 
correspond to the multipoles I — 1952 and t = 2635, respectively. The contours from the copula likelihood are nearly 
identical to those from the Gaussian likelihood, because the % 2 distribution approaches to Gaussian distribution at 
these small scales due to the central limit theorem. We can see a positive strong correlation which is mainly attributed 
to the fact that the nonlinear gravitational evolution causes the non-Gaussian contribution and correlations between 
different multipole bins (e.g., see in Fig. 7 in (l8j|). 

To confirm a much better reliability of our Copula likelihood, we examine how the Copula likelihood is statistically 
better relative to the Gaussian likelihood based on the information criterion theory. We use the Akaike Information 
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Criterion (AIC) [2!| [3^, to evaluate which model is more preferable. The AIC is defined as 

AIC = -21n£ max + 2fc, (31) 

where £ max is the maximum likelihood achievable by the model under a certain dataset, and k is the number of 
free parameters of the model. In the present analysis £ max is directly obtained by substituting P K (£i) into Eqs. (|25j) 
and (|26p at each realization, because model parameters have been already fixed to give the maximum probability for 
the 1000 realized simulation data. The meaning of the AIC is clearly understood as an extension of the maximum 
likelihood method. An explanation for astronomers can be found in |29| and there are also useful applications of the 
AIC in their paper. Using the AIC enables us to compare the goodness of a certain model with that of another type 
directly. Figure @] shows the PDF of the AIC difference, AAIC = AlCcaussian — AlCcopuia among 1000 realizations. 
We see that the AIC is significantly reduced by the Copula likelihood (positive AAIC) for the vast majority of the 
realizations. This means that the Copula likelihood function represents the distribution of the convergence power 
spectra from Monte Carlo simulations better than the Gaussian likelihood one. We estimate that the mean value 
of the AIC difference is 2.35. Some rough rules of thumb are available and useful for estimation of the goodness 
of models. As shown in [35|, the mean value of the AIC difference we have estimated is regarded as Considerably 
less which means the evidence to support the Gaussian likelihood is considerably less than the Copula likelihood. 
Therefore we can conclude that our Copula likelihood is strongly favored compared to the Gaussian likelihood for 
cosmic shear power spectrum. 

Finally, we discuss how the difference between the two likelihoods affects the cosmological parameter estimation, 
paying attention to the number of bins. Since weak lensing power spectra are expected to be very smooth in £, it is very 
likely that little information is lost by binning as long as the bins are narrow compared to the width of any features. 
If we consider an observation over the sky coverage / s k y (/sky = f2 s /47r) with useful signal at if ^ I ;$ £ max , the total 
number of modes one can obtain is estimated to be n t ~ /sky(^ ma x — ^/)- If the number of bins is small enough, 
the distribution of the power at each bin can be approximated by a Gaussian because a sufficient number of modes 
are in that bin. As is shown in [36] if the number of bins satisfies a condition Uf, <C / s Yy ^max assuming ^ ax S> l 2 j , 
the Gaussian approximation for the one-point PDF at each bin becomes very good. In that case one can use the 
multivariate Gaussian likelihood in order to estimate the cosmological parameters, which will simplify the parameter 
estimation analysis. Now let us apply this argument to our examples. In our example with / s k y = 6.25 x 10~ 4 and 

1/2 

f max ~ 1000, the square root of the total number of the modes is n t ~ / sky dax = 25. This number is not so large 
compared to the number of bins, i.e., n& — 13. Therefore, in this case an accurate likelihood function should be 
used in order toget unbiased cosmological parameter constraints, instead of the conventional multivariate Gaussian 
likelihood [see, |26j. Meanwhile, if we consider a future type survey which has / s k y = 0.05 and £ max ~ 1000 with the 

1/2 

same number of bins, the square root of the total number of the modes is / sky ^ max = 224 which is much larger than 
the number of bins. Then, the Gaussian approximation could be fine in this case [2(| 5 . Note, however, that Gaussian 



This is a rough discussion because our binning is logarithmic in I space and each bin does not contain an equal number of the modes. 
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FIG. 4: Probability distribution of the Akaike Information Criterion (AIC) difference between the Copula likelihood (Eq. 
and the Gaussian likelihood (Eg. I26[) . among 1000 realizations. 



approximation may be violated by natural binning which has optimal £ resolution, regardless of how £ max value is. 
Therefore, we suggest the Copula likelihood should be used when an optimal binning is done to keep the cosmological 
information as much as possible. The impact of the difference between the two likelihoods on cosmological parameter 
estimations is illustrated and discussed carefully in a companion paper [26| . 



V. CONCLUSION 



It has been becoming an important issue to obtain not only accurate statistics such as a power spectrum and 
covariance matrix but also an accurate likelihood function for the precision cosmology. In this work, we introduced 
a statistical tool called a copula into cosmology in a rather pedagogical way. The copula is a function to generate 
an n-point CDF from the given one-point CDFs and prescribed dependence structure of variables. We then applied 
the copula to the cosmic convergence power spectrum estimated from 1000 realizations obtained from ray-tracing 
simulations generated by (l8| . 

By taking into account the fact that the one-point PDF of the convergence power spectrum is well approximated by 
\ 2 distribution, we showed that the Gaussian copula can reproduce the n-dimensional shape of the likelihood function 
of the convergence power spectrum better than the multivariate Gaussian likelihood, for the simulated 13-binned data 
expected from the assumed survey area of i7 s = 25 deg 2 . This is a main result of the present paper. The differences 
between the two likelihood models become significant at lowest multipole bins. 

The deviation from the multivariate Gaussian will heavily depend on the width of the binning. We discussed that 
a sufficiently sparse binning will make it possible to use a multivariate Gaussian likelihood for the future cosmic shear 
survey. However, if one makes an optimal binning (A£ ~ 1) to keep the information as much as possible, non-Gaussian 
corrections will become important. In this case the copula will provide us an appropriate likelihood function in a 
convenient way by relating the one-point CDFs to the n-point CDF. 
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